Fully selfconsistent GW calculations for molecules 
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We calculate single-particle excitation energies for a series of 33 molecules using fully selfcon- 
sistent GW, one-shot Go Wo, Hartree-Fock (HF), and hybrid density functional theory (DFT). All 
calculations are performed within the projector augmented wave (PAW) method using a basis set 
of Wannier functions augmented by numerical atomic orbitals. The GW self-energy is calculated 
on the real frequency axis including its full frequency dependence and off-diagonal matrix elements. 
The mean absolute error of the ionization potential (IP) with respect to experiment is found to be 
4.4, 2.6, 0.8, 0.4, and 0.5 eV for DFT-PBE, DFT-PBEO, HF, G W [HF], and selfconsistent GW, 
respectively. This shows that although electronic screening is weak in molecular systems its inclu- 
sion at the GW level reduces the error in the IP by up to 50% relative to unscreened HF. In general 
GW overscreens the HF energies leading to underestimation of the IPs. The best IPs are obtained 
from one-shot Go Wo calculations based on HF since this reduces the overscreening. Finally, we find 
that the inclusion of core- valence exchange is important and can affect the excitation energies by as 
much as 1 eV. 

PACS numbers: 31.15.A-,33.15.Ry,31.15.V- 



I. INTRODUCTION 

Density functional theory (DFT)i with the single- 
particle Kohn-Sham (KS) scheme^ is today the most 
widely used approach to the electronic structure problem 
of real materials in both solid state physics and quantum 
chemistry. While properties derived from total energies 
(or rather total energy differences) are accurately pre- 
dicted by DFT, it is well known that DFT suffers from 
a band gap problem implying that the single-particle KS 
eigenvalues cannot in general be interpreted as real quasi- 
particle (QP) excitation energies. In particular, semilocal 
exchange-correlation functionals severely underestimate 
the fundamental gap of both insulators, semi-conductors, 
and molecules 

The hybrictfr— and screened hybrid^ functionals, 
which admix around 25% of the (screened) Fock ex- 
change with the local DFT exchange, generally improve 
the description of band gaps in bulk semi-conductors 
and insulators^. However, the orbital energies obtained 
for finite systems using such functionals still underes- 
timate the fundamental gap, I p — E a , (the difference 
between ionization potential and electron affinity) by 
up to several electron volts. In fact, for molecules the 
pure Hartree-Fock (HF) eigenvalues are usually closer 
to the true electron addition/removal energies than are 
the hybrid DFT eigenvalues. This is because HF is 
self-interaction free and because screening of the ex- 
change interaction is a relatively weak effect in molecular 
systems ^ 11 ' 13 . On the other hand, in extended systems 
the effect of self-interaction is less important and the 
long range Coulomb interaction becomes short ranged 
due to dynamical screening. As a consequence HF breaks 
down in extended systems leading to dramatically over- 
estimated band gaps and a qualitatively incorret descrip- 
tion of metalsJ^— 



The many-body GW approximation of Hedinil 
has been widely and succesfully used to calculate 
QP band structures in metals, semi-conductors, and 
insulators jiL^r— The GW approximation can be viewed 
as HF with a dynamically screened Coulomb interac- 
tion. The fact that the screening is determined by the 
system itself instead of being fixed a priori as in the 
screened hybrid schemes, suggests that the GW method 
should be applicable to a broad class of systems rang- 
ing from metals with strong screening to molecules with 
weak screening. With the entry of nanoscience the use of 
GW has been extended to low-dimensional systems and 
nanostructures^ir— and more recently even nonequilib- 
rium phenomena like quantum transporter—. In view 
of this trend it is important to establish the perfor- 
mance of the GW approximation for other systems than 
the crystalline solids. In this work we present first- 
principles benchmark GW calculations for a series of 
small molecules. In a closely related study we compared 
GW and Hartree-Fock to exact diagonalization results for 
semi-empirical PPP models of conjugated molecules 3 ^. 
The main conclusions from the two studies regarding the 
qualities of the GW approximation in molecular systems, 
are very consistent. 

Most GW calculations to date rely on one or several 
approximations of more technical character. These in- 
clude the plasmon pole approximation, the linearized 
QP equation, neglect of off-diagonal matrix elements 
in the GW self-energy, analytic continuations from the 
imaginary to the real frequency axis, neglect of core 
states contributions to the self-energy, neglect of self- 
consistency. The range of validity of these approxima- 
tions has been explored for solid state systems by a num- 
ber of authors^—, however, much less is known about 
their applicability to molecular systems^. Our imple- 
mentation of the GW method avoids all of these tech- 
nical approximations allowing for a direct and unbiased 
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assessment of the GW approximation itself. 

Here we report on single-shot Go Wo and fully self- 
consistent GW calculations of QP energies for a set of 
33 molecules. The calculated IPs are compared with ex- 
perimental values as well as single-particle eigenvalues 
obtained from Hartree-Fock and DFT-PBE /PBEO theo- 
ries. As additional benchmarks we compare to second- 
order Moller-Plesset (MP2) and DFT-PBE total energy 
differences between the neutral and cation species. Spe- 
cial attention is paid to the effect of sclfconsistency in the 
GW self-energy and the role of the initial Green func- 
tion, Go, used in one-shot Go Wo calculations. The use 
of PAW rather than pseudopotentials facilitate the in- 
clusion of core- valence exchange, which we find can con- 
tribute significantly to the HF and GW energies. Our 
results show that the GW approximation yields accurate 
single-particle excitation energies for small molecules im- 
proving both hybrid DFT and full Hartree-Fock results. 

The paper is organized as follows. In Sec. |H] we de- 
scribe the theoretical and numerical details behind the 
GW calculations, including the augmented Wannier func- 
tion basis set, the self-consistent solution of the Dyson 
equation, and the evaluation of valence-core exchange 
within PAW. In Sec. IIIII we discuss and compare the 
results of G W , GW, HF, PBEO, and PBE calculations. 
We analyze the role of dynamical screening, and discuss 
the effect of self-consistency in the GW self-energy. We 
conclude in Sec. [IVJ 

II. METHOD 

A. Augmented Wannier function basis 

For the GW calculations we apply a basis set consist- 
ing of projected Wannier functions (PWF) augmented 
by numerical atomic orbitals (NAO). The PWFs, 0i, are 
obtained by maximizing their projections onto a set of 
target NAOs, &Aim, subject to the condition that they 
span the set of occupied cigenstates, ip n . Thus we maxi- 
mize the functional 

^ = E E i(^i^ m )i 2 a) 

i A,l,m 

subject to the condition span{^} D span{f/; Tl } occ as 
described in Ref. |U The target NAOs are given by 
®Ai m (r) = (Ai(r)Y hn (r) where ( A i is a modified Gaus- 
sian which vanish outside a specified cut-off radius, and 
Yim are the spherical harmonics corresponding to the va- 
lence of atom A. The number of PWFs equals the num- 
ber of target NAOs. For example we obtain one PWF 
for H (/ max = 0), and four PWFs for C (Z max = 1). The 
PWFs mimick the target atomic orbitals but in addition 
they allow for an exact representation of all the occupied 
molecular eigenstates. The latter are obtained from an 
accurate real-space PAW-PBE calculatio n 4 ^ 46 . 

The PWFs obtained in this way provide an exact rep- 
resentation of the occupied PBE eigenstates. However, 



this does not suffice for GW calculations because the po- 
larizability, P, and the screened interaction, W, do not 
live in this subspace. Hence we augment the PWFs by 
additional NAOs including so-called polarization func- 
tions which have / = l max + 1 and/or extra radial func- 
tions (zeta functions) for the valence atomic orbitals. For 
more details on the definition of polarization- and higher 
zeta functions we refer to Ref. |46| . To give an exam- 
ple, a double-zeta-polarized (DZP) basis consists of the 
PWFs augmented by one set of NAOs corresponding to 
I = 0, l m &x and one set of polarization orbitals. Note 
that the notation, SZ, SZP, DZ, DZP, etc., is normally 
used for pure NAO basis sets, but here we use it to de- 
note our augmented Wannier basis set. We find that the 
augmented Wannier basis is significantly better for HF 
and GW calculations than the corresponding pure NAO 
basis. 

The GW and HF calculations presented in Sec. IIIII 
were performed using a DZP augmented Wannier basis. 
This gives a total of 5 basis functions per H, Li, and Na, 
and 13 basis functions for all other chemical elements con- 
sidered. In Sec. IIII CI we discuss convergence of the GW 
calculations with respect to the size of the augmented 
Wannier basis. 



B. GW calculations 

The HF and GW calculations for isolated molecules 
are performed using a Green function code developed 
for quantum transport^ In principle, this scheme is de- 
signed for a molecule connected to two electrodes with 
different chemical potentials /i£ and fj,R. However, the 
case of an isolated molecule can be treated as a special 
case by setting /i£ = fiR = /i and modelling the cou- 
pling to electrodes by a constant imaginary self-energy, 
^l/r = iV' The chemical potential /i is chosen to lie in 
the HOMO-LUMO gap of the molecule and the size of 
77, which provides an artifical broadening of the discrete 
levels, is reduced until the results have converged. In this 
limit of small r\ the result of the GW calculation becomes 
independent of the precise position of /i inside the gap. 

In Ref. H3 the GW-transport scheme was described for 
the case of an orthogonal basis set and for a truncated, 
two-index Coulomb interaction. Below we generalize the 
relevant equations to the case of a non-orthogonal basis 
and a full four-index Coulomb interaction. Some cen- 
tral results of many-body perturbation theory in a non- 
orthogonal basis can be found in Ref. 

The central object is the retarded Green function, G r , 

G r (e) = {(e + i^S-IlKs + v^-AvK-KclG^T 1 (2) 

In this equation all quantities are matrices in the aug- 
mented Wannier basis, e.g. H^,ij = {4>i\H-Ks\4>j) is 
the KS Hamiltonian matrix and SV, — ((f>i\(f>j} is an 
overlap matrix. The term A«h represents the change 
in the Hartree potential relative to the DFT Hartree 
potential already contained in -ffKSi see Appendix [X] 
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The local xc-potential, u xc , is subtracted to avoid dou- 
ble counting when adding the many-body self-energy, 
X XC [G]. As indicated, the latter depends on the Green 
function and therefore Eq. ([5} must in principle be solved 
self-consistently in conjuction with the self-energy. 

In the present study S xc is either the bare exchange 
potential or the GW self-energy. To be consistent with 
the code used for the calculations, we present the equa- 
tions for the GW self-energy on the so-called Keldysh 
contour. However, under the equilibrium conditions con- 
sidered here the Keldysh formalism is equivalent to the 
ordinary time-ordered formalism. 

The GW self-energy is defined by 

Sg w (r,r') = iJ2Gkl(ry + )W ikdl (T,r'), (3) 

kl 

where r and r' are times on the Keldysh contour, C. 
The dynamically screened Coulomb interaction obeys the 
Dyson-like equation 

W ijM (T,T') = V ijM 6c(T,T') 

,pqPpq,rs("T, T\ ) W ra ^l {j\ , t') , (4) 

pqrs C 

and the polarization bubble is given by 

fy,w(r,r') = -iG tk (T,T')G l3 (T',T). (5) 

In the limit of vanishing polarization, P = 0, W reduces 
to the bare Coulomb interaction 

= // |^f^W^W«(0^(0 (6) 

and the GW self-energy reduces to the exchange potential 
of HF theory. 

From the above equations for the contour-ordered 
quantities, the corresponding real time components, i.e. 
the retarded, advanced, lesser, and greater components, 
can be obtained from standard conversion rule a 49 i 50 . For 
completeness we give the expressions for the real time 
components of the GW equations in Appendix [SJ 

The time/energy dependence of the dynamical quan- 
tities G, W, P, and E, is represented on a uniform grid. 
We switch between time and energy domains using the 
Fast Fourier Transform in order to avoid time consuming 
convolutions. A typical energy grid used for the GW cal- 
culations in this work ranges from -150 to 150 eV with 
a grid spacing of 0.02 eV. The code is parallelized over 
basis functions and energy grid points. We use a Pulay 
mixing scheme for updating the Green function G r when 
iterating Eq. to self-consistency as described in Ref. 

S3. 

We stress that no approximation apart from the fi- 
nite basis set is made in our implementation of the GW 
approximation. In particular the frequency dependence 
is treated exactly and analytic continuations from the 
imaginary axis are avoided since we work directly on the 
real frequency/time axis. The price we pay for this is the 
large size of the energy grid. 



C. Spectral function 

The single-particle excitation spectrum is contained in 
the spectral function 

A(e) = l (G r (e)-[G r (e)}^. (7) 

For a molecule A(e) shows peaks at the QP energies 
e n = E n (N + 1) - E (N) and e n = E (N) - E n (N - 1) 
corresponding to electron addition and removal energies, 
respectively. Here E n (N) denotes the energy of the nth 
excited state of the system with N electrons and N refers 
to the neutral state. 

When the Green function is evaluated in a non- 
orthogonal basis, like the augmented Wannier basis used 
here, the projected density of states for orbital & be- 
comes 

Di(e) = [SA(e)S] u / 2irS u , (8) 

where matrix multiplication is implied^ Correspond- 
ingly, the total density of states, or quasiparticle spec- 
trum, is given by 

D(e) = Tv(A(e)S)/2n. (9) 

D. Calculating Coulomb matrix elements 

The calculation of all of the Coulomb matrix elements, 
Vij.kh is prohibitively costly for larger basis sets. For- 
tunately the matrix is to a large degree dominated by 
negligible elements. To systematically define the most 
significant Coulomb elements, we use the product basis 
technique of Aryasetiawan and Gunnarsso n 51 ' 52 . In this 
approach, the pair orbital overlap matrix 

Sij,ki = (riij\n k i), (10) 

where riy(r) = <fi*(r)<f)j(r) is used to screen for the sig- 
nificant elements of V . 

The eigenvectors of the overlap matrix Eq. (fIT)|) rep- 
resents a set of "optimized pair orbitals" and the eigen- 
values their norm. Optimized pair orbitals with insignif- 
icant norm must also yield a reduced contribution to the 
Coulomb matrix, and are omitted in the calculation of 
V. We limit the basis for V to optimized pair orbitals 
with a norm larger than 10 _5 q,q~ 3 . This gives a signifi- 
cant reduction in the number of Coulomb elements that 
needs to be evaluated, and it reduces the matrix size of 
P(e) and W(e) correspondingly, see Appendix [Al 

The evaluation of the double integral in Eq. © is 
efficiently performed in real space by solving a Poisson 
equation using multigrid technique s 45 ' 53 . 

E. Valence-core exchange 

All inputs to the GW/HF calculations, i.e. the selfcon- 
sistent Kohn-Sham Hamiltonian, -HkSj the xc potential 
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v xc , the Coulomb matrix elements, Vij t ki, are calculated 
using the real-space PAW-2 4 . code GPAW— i 4 ^. 

In GPAW, the core electrons (which are treated scalar- 
relativistically) are frozen into the orbitals of the free 
atoms, and the Kohn-Sham equations are solved for the 
valence states only. Unlike pseudo potential schemes, 
these valence states are subject to the full potential of 
the nuclei and core electrons. This is achieved by a parti- 
tioning scheme, where quantities are divided into pseudo 
components augmented by atomic corrections. The oper- 
ators obtained from GPAW are thus full-potential quan- 
tities, and the wave functions from which the Wannier 
basis functions are constructed correspond to the all- 
electron valence states. Ref. HH describes how the all- 
electron Coulomb elements can be determined within the 
PAW formalism. 

Since both core and all-electron valence states are 
available in the PAW method, we can evaluate the contri- 
bution to the valence exchange self-energy coming from 
the core electrons. As the density matrix is simply the 
identity matrix in the subspace of atomic core states, this 
valence-core exchange reads 



y^corc 
x.ij 



core 

£ 

k 



yikjjk ? 



(ii) 



where i,j represent valence basis functions. We limit 
the inclusion of valence-core interactions to the exchange 
potential, neglecting it in the correlation. This is reason- 
able, because the polarization bubble, P, involving core 
and valence states will be small due to the large energy 
difference and small spatial overlap of the valence and 
core states. This procedure was used and validated for 
solids in Ref. M We find that the elements of £™J| 
can be significant - on average 1.2 eV for the HOMO - 
and are larger (more negative) for the more bound or- 
bitals which have larger overlap with the core states. In 
general, the effect on the HOMO-LUMO gap is to en- 
large it, on average by 0.4 eV because the more bound 
HOMO level is pushed further down than the less bound 
LUMO state. In the case of solids, the role of valence- 
core interaction has been investigated by a number of 
authors ^ 42 ' 55 . Here the effect on the QP band gap 
seems to be smaller than what we find for the molec- 
ular gaps. We note that most GW calculations rely on 
pseudopotential schemes where these valence-core inter- 
actions are not accessible. In such codes, the xc contri- 
bution from the core electrons are sometimes estimated 
by £ X ° IC ~ ^xcM — w xc [n V ai] where rt V ai is the valence 
electron density, but as the local xc potential is a non- 
linear functional of the density, this procedure is not well 
justified. Instead we subtract the xc potential of the full 
electon density n, and add explicitly the exact exchange 
core contribution. 



III. RESULTS 

In Fig. [T] we compare the calculated HOMO en- 
ergies with experimental ionization potentials for the 
33 molecules listed in Table HI The geometries of the 
molecules, which all belong to the G2 test set, are taken 
from Ref. [56|. The different HOMO energies correspond 
to: DFT-PBE^ and DFT-PBEOl eigenvalues, Hartree- 
Fock eigenvalues, and fully selfconsistent GW. The GW 
energies are obtained from the peaks in the correspond- 
ing density of states Eq. ^ extrapolated to r\ = (rj 
gives an artificial broadening of the delta peaks). 
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FIG. 1: (Color online) Calculated negative HOMO energy 
versus experimental ionization potential. Both PBE and 
PBEO systematically understimates the ionization energy due 
to self-interaction errors while HF overestimates it slightly. 
The dynamical screening from the GW correlation lowers the 
HF energies bringing them closer to the experimental values. 
Numerical values are listed in Table [I] 

We stress the different meaning of fully selfconsistent 
GW and the recently introduced method of quasiparti- 
cle selfconsistent GW— : In fully selfconsistent GW the 
Green function obtained from Dyson's equation Eq. ([2]) 
with £;e C [G] = Egw[G] is used to calculate the £gw of 
the next iteration. In QP-selfconsistent GW, Sqw is al- 
ways evaluated using a non-interacting Green function 
and the self-consistency is obtained when the difference 
between the non- interacting GF and the interacting GF, 
is minimal. 

Fig. Q] clearly shows that both the PBE and PBEO 
eigenvalues of the HOMO severely underestimates the 
ionization potential. The average deviation from the ex- 
perimental values are 4.35 eV and 2.55 eV, respectively. 
The overestimation of the single-particle eigenvalues of 
occupied states is a well known problem of DFT and can 
be ascribed to the insufficient cancellation of the self- 
interaction in the Hartree potential! 4 ^ Part of this self- 
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interaction is removed in PBEO. However, the fact that 
the HF results are significantly closer to experiments in- 
dicates that the 25% Fock exchange included in the PBEO 
is not sufficient to cure the erroneous description of (oc- 
cupied) molecular orbitals. On the other hand PBEO 
gives good results for band gaps in semi-conductors and 
insulators where in contrast full Hartree-Fock does not 
perform wellJ^r— We conclude that the amount of Fock 
exchange to be used in the hybrid functionals to achieve 
good quasiparticle energies is highly system dependent. 
A similar problem is encountered with self-interaction 
corrected exchange-correlation functionals. 13 




Ae HF (eV) 

FIG. 2: The deviation of the calculated HOMO energy from 
the experimental ionization potential in GW and HF, respec- 
tively. The vertical displacement of points from the line x — y 
gives the difference between the GW and HF energies and rep- 
resents the effect of screening. Notice that the GW correction 
is always negative (corresponding to higher HOMO energy) 
and that it generally overcorrects the HF energies. Also no- 
tice that the GW correction is larger for molecules where HF 
presents the largest overestimation of the ionization potential. 

As can be seen from Fig. [T] GW performs better than 
Hartree-Fock for the HOMO energy yielding a mean ab- 
solute error with respect to experiments of 0.5 eV com- 
pared to 0.81 eV with Hartree-Fock. As expected the dif- 
ference between HF and GW is not large on an absolute 
scale (around 1 eV on average, see Table [Hi illustrating 
the fact that screening is weak in small molecules. On a 
relative scale selfconsistent GW improves the agreement 
with experiments by almost 30% as compared to HF. 

To gain more insight into the influence of screening on 
the orbital energies, we compare in Fig. [2] the deviation 
of the HF and GW energies from IP cxp . The GW self- 
energy can be split into the bare exchange potential and 
an energy-dependent correlation part 

E GW (r, r'; e) = v x {r, r') + S corr (r, r'; e) (12) 



Accordingly the quasiparticle energy can be written as 
the bare HF energy and a correction due to the energy- 
dependent part of the GW self-energy (the dynamical 
screening term) 

sT = e™ + A« w (13) 

In Fig. [2] the line y — x corresponds to A^ w = 0, and 
the vertical displacement from the line thus represents 
the effect of screening on the calculated HOMO energy. 
We first notice that the effect of screening is to shift the 
HOMO level upwards in energy, i.e. to reduce the ion- 
ization potential. This can be understood by recalling 
that the Hartree-Fock eigenvalue represents the energy 
cost of removing an electron from the HOMO when or- 
bital relaxations in the final state are neglected (Koop- 
mans' theorem^). In Ref. 37] we showed, on the basis 
of GW and exact calculations for semi-empirical models 
of conjugated molecules, that A„ w mainly describes the 
orbital relaxations in the final state and to a lesser extent 
accounts for the correlation energy of the initial and final 
states. This explains the negative sign of A^ w because 
the inclusion of orbital relaxation in the final state lowers 
the energy cost of removing an electron. We note that 
this is different from the situation in extended, periodic 
systems where orbital relaxations vanish and the main ef- 
fect of the GW self-energy is to account for correlations 
in the initial and final states. 

In TableUwe list the calculated HOMO energy for each 
of the 33 molecules. In addition to selfconsistent GW we 
have performed one-shot Go Wo calculations based on the 
HF and PBE Green's function, respectively. The best 
agreement with experiment is obtained for GoWo[HF]. 
This is because the relatively large Hartree-Fock HOMO- 
LUMO gap reduces the (over-)screening described by the 
resulting GW self-energy. There are not many GW cal- 
culations for molecules available in the litterature. Below 
Table U we list the few we have found. As can be seen 
they all compare quite well with our results given the 
differences in the implementation of the GW approxima- 
tion. 

For comparison we have included the HOMO energy 
predicted by second order M0ller-Plesset theory (MP2) 
[taken from Ref. [5!| with a Gaussian 6-311G** basis 
set. These are generally very close to our calculated HF 
values, with a tendency to lower energies which worsens 
the agreement with experiment slightly as compared to 
HF. 

We have also calculated the DFT-PBE total energy dif- 
ference between the neutral and cation species, E(N) — 
E(N — 1), see last column of Table |TJ This procedure 
leads to IPs in very good agreement with the experimen- 
tal values (MAE of 0.24 eV). We stress that although 
this method is superior to the GW method for the IP of 
the small molecules studied here, it can yield only the 
HOMO and LUMO levels while higher excited states are 
inaccessible. Moreover it applies only to isolated systems 
and cannot be directly used to probe QP levels of e.g. a 
molecule on a surface. 
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TABLE I: Experimental ionization potential (first column) and HOMO energy calculated using different approximations for 
exchange and correlation. "X-eig" refers to a single-particle eigenvalue while "X-tot" refers to a total energy difference, 
E(N) — E(N — 1). The GoWo(PBE) energies have been obtained from the QP equation while the GW and Go Wo energies are 
obtained from the DOS in Eq. {5). Last row shows the mean absolute error (MAE) with respect to experiments. All energies 
are in eV. 



Molecule 


Expt. (o) 


PBE-eig 


PBEO-eig 


HF-eig 


GW 


GoWo(HF) 


G W (PBE)-QP 


MP2 (a) 


PBE-tot 


LiH 


7.90 


4.34 


5.81 


8.14 


8.0 (i,) 


8.2< 6 > 


8.0 


8.20 


8.02 


Li2 


5.11 


2.96 


3.62 


4.62 


4.6 


4.7 


4.4 


4.91 


5.09 


LiF 


11.30 


6.00 


8.62 


13.26 


11.7 


11.2 


12.0 


12.64 


11.87 


Na 2 


4.89 


2.81 


3.38 


4.16 


4.1 


4.3 


4.7 


4.48 


4.97 


NaCl 


9.80 


5.24 


6.92 


9.78 


9.0 


9.2 


8.8 


9.63 


9.37 


CO 


14.01 


9.05 


10.98 


14.80 


13.4 


14.1 


13.9 


15.08 


13.88 


COa 


13.78 


9.08 


11.09 


14.50 


13.1 


13.3 


13.6 


14.71 


13.64 


CS 


11.33 


7.40 


9.09 


12.31 


10.8 


11.7 


11.0 


12.58 


11.31 


C2H2 


11.49 


7.20 


8.64 


11.05 


10.6 


11.1 


11.2 


11.04 


11.39 


C2H4 


10.68 


6.79 


8.11 


10.11 


9.8 


10.4 


9.6 


10.18 


10.67 


CH4 


13.60 


9.43 


11.29 


14.77 


14.1 


14.4 


14.4 (c) 


14.82 


14.10 


CH3CI 


11.29 


7.08 


8.80 


11.68 


11.0 


11.4 


11.1 


11.90 


11.10 


CH3OH 


10.96 


6.31 


8.49 


12.14 


10.7 


10.8 


10.5 


12.16 


10.72 


CH 3 SH 


9.44 


5.60 


7.09 


9.50 


8.8 


9.0 


8.4 


9.73 


9.29 


Cl 2 


11.49 


7.32 


9.02 


12.03 


10.9 


11.3 


11.5 


12.37 


11.22 


C1F 


12.77 


7.90 


9.88 


13.33 


12.4 


12.4 


13.0 


13.63 


12.48 


F 2 


15.70 


9.43 


12.42 


17.90 


15.2 


15.2 


16.2 


18.20 


15.39 


HOC1 


11.12 


6.68 


8.66 


11.93 


10.6 


10.8 


11.0 


12.23 


10.95 


HC1 


12.74 


8.02 


9.78 


12.96 


12.2 


12.5 


12.5 


13.02 


12.71 


H2O2 


11.70 


6.38 


8.78 


13.06 


11.0 


11.1 


11.1 


13.00 


11.18 


H 2 CO 


10.88 


6.28 


8.37 


11.93 


10.4 


10.5 


10.6 


11.97 


10.80 


HCN 


13.61 


9.05 


10.67 


13.19 


12.7 


13.2 


12.4 


13.33 


13.67 


HF 


16.12 


9.61 


12.47 


17.74 


16.0 


15.6 


15.7 


17.35 


16.27 


H 2 


12.62 


7.24 


9.59 


13.88 


12.3 


12.1 


11.9 (d) 


13.62 


12.88 


NH 3 


10.82 


6.16 


8.11 


11.80 


10.8 


11.0 


10.6 


11.57 


11.02 


N 2 


15.58 


10.28 


12.51 


16.21 


15.1 


15.7 


15.6 


16.41 


15.39 


N 2 H 4 


8.98 


5.75 


7.67 


11.06 


9.8 


10.1 


9.5 


11.07 


9.90 


SH 2 


10.50 


6.29 


7.79 


10.48 


9.8 


10.1 


9.9 


10.48 


10.38 


SO2 


12.50 


8.08 


9.96 


13.02 


11.3 


11.7 


11.7 


13.46 


12.12 


PH 3 


10.95 


6.79 


8.17 


10.38 


9.9 


10.3 


10.0 


10.50 


10.39 


P2 


10.62 


7.09 


8.21 


9.65 


9.2 


9.8 


9.0 


10.09 


10.37 


SiH 4 


12.30 


8.50 


10.13 


12.93 


12.3 


12.6 


12.4 (e) 


13.25 


11.95 


Si 2 H 6 


10.53 


7.27 


8.54 


10.82 


10.2 


10.6 


9.9 


11.03 


10.36 


SiO 


11.49 


7.46 


9.14 


11.78 


10.9 


11.2 


11.3 


11.82 


11.27 


MAE 




4.35 


2.55 


0.81 


0.5 


0.4 


0.5 


0.82 


0.24 



^From Ref. IBISI . The MP2 calculations use a Gaussian 6-311G** basis set. 

(b) To be compared with the GW value 7.85 and the G W (HF) value 8.19 reported in Ref. H 

(c) To be compared with the G W (LDA) value 14.3 reported in Ref. Hi]. 

(d) To be compared with the G W (LDA) value 11.94 reported in Ref. [H 

^ e 'To be compared with the GoWo(LDA) values 12.7 and 12.66 reported in Refs. [21] andl22l. respectively. 



In Table |TT] we provide an overview of the compara- are the same as those listed in the last row of Table HH 
tive performance of the different methods. Shown is the 
mean average deviation between the IPs calculated with 
the different methods as well as the experimental values. 
Note that the numbers in the experiment row/column 
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TABLE II: Mean absolute deviation between the IPs of the 33 molecules calculated with the different methods and experiment. 
The mean absolute deviation with respect to experiment coincide with the last row in Table [J 



Method 


Expt. (a) 


PBE-eig 


PBEO-eig 


HF-eig 


GW 


G Wo[HF] 


MP2 (a) 


PBE-tot 




u.uu 


A 


A. DO 


n si 

U.ol 


u.o 




n so 




PBE 


4.35 


0.00 


1.79 


4.90 


3.9 


4.1 


4.99 


4.27 


PBEO 


2.55 


1.79 


0.00 


3.11 


2.1 


2.3 


3.20 


2.48 


HF 


0.81 


4.90 


3.11 


0.00 


1.0 


0.8 


0.17 


0.80 


GW 


0.5 


3.9 


2.1 


1.0 


0.00 


0.3 


1.1 


0.4 


G Wo[HF] 


0.4 


4.1 


2.3 


0.8 


0.3 


0.00 


0.9 


0.3 


MP 2 


0.82 


4.99 


3.20 


0.17 


1.1 


0.9 


0.00 


0.84 


PBE-tot 


0.24 


4.27 


2.48 


0.80 


0.4 


0.3 


0.84 


0.00 



(a) Data taken from Ref. H 
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FIG. 3: (Color online) Density of states for the NH3 molecule 
calculated in HF and GW, respectively. Arrows mark the level 
corresponding to the HOMO in the two calculations. The in- 
tersection between the line y = e — e™ and the real part of 
(^HOMol s co rr (e)|i/'HOMo) (green curve) determines the posi- 
tion of the GW level. 



A. Linearized quasiparticle equation 

In the conventional GW method the full Green func- 
tion of Eq. (J2) is not calculated. Rather one obtains the 
quasiparticle energies from the quasiparticle equation 

£ Q p = e° n + Z„«|£ GW (4) - v^° n ). (14) 

where ip^ and e„ are eigenstates and eigenvalues of an 
approximate single-particle Hamiltonian (often the LDA 
Hamiltonian), and 



DFT I 



ds 



(15) 



Moreover the GW self-energy is evaluated non- 
sclfconsistently from the single-particle Green function, 
i.e. Sqw = iG W[G ], with G (z) = {z- Ho)- 1 . 

The quasiparticle equation (I14[) relies on the assump- 
tion that off-diagonal matrix elements, (ip®\T,Gw(£n) ~ 



fxcIV"™); can be neglected, and that the frequency de- 
pendence of Sqw can be approximated by its first or- 
der Taylor expansion in a sufficiently large neighbor- 
hood of e". We have found that these two assumptions 
are indeed fullfilled for the molecular systems studied 
here. More precisely, for the GW and G W (HF) self- 
energies, the QP energies obtained from Eq. (fT4"|) are 
always very close to the peaks in the density of states 
Eq. ([9]). We emphasize that this result could well be re- 
lated to the rather large level spacing of small molecules, 
and may not hold for extended systems. An example is 
presented in Fig. [3] which shows the full HF and GW 
density of states for NH3 together with the real part of 



HOMOl^ corr 



HOMO/ 



As explained in the follow- 
ing section this is not quite the case for the Go Wo (PBE) 
calculations. 



B. Go-dependence 

As stated in the previous section the GW and 
GoWo(HF) energies can be obtained cither from the full 
spectral function or from the QP equation. In this case, 
returning to Table Q] we see that Go Wo (HF) yields sys- 
tematically larger IPs than GW. This is easy to under- 
stand since Ghf describes a larger HOMO-LUMO gap 
than Ggwj and therefore produces less screening. When 
the PBE rather than the HF Green function is used to 
evaluate the GW self-energy, we find that the spectral 
function obtained from Eq. @ does not resemble a 
simple discrete spectrum. In fact the peaks are signif- 
icantly broadened by the imaginary part of Eqw and it 
becomes difficult to assign precise values to the QP ener- 
gies. Apart from the spectral broadening, the molecular 
gap is significantly reduced with respect to its value in the 
GW and GoWo(HF) calculations. Both of these effects 
are due to the very small HOMO-LUMO gap described 
by Gpbe which leads to severe overscreening and QP life- 
time reductions. A similar effect was observed by Ku and 
Eguiluz in their comparison of GW and Go Wo (LDA) for 
Si and Ge crystals^. 

The problems encountered when attempting to solve 
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the Dyson equation ([2]) using the Go Wo (PBE) self-energy 
occur due to the large mismatch between e^ BE and e^ p . 
On the other hand, in the QP equation, the GW self- 
energy is evaluated at e° rather than e^ p . As a conse- 
quence the unphysical broadening and overscreening is 
avoided and a well defined QP energy can be obtained 
(last column in Table [I]). 

To summarize, Go can have a very large effect on the 
QP spectrum when the latter is obtained via the Dyson 
equation ([2]). In particular, the use of a Go with a too 
narrow energy gap (as e.g. the Gpbe) can lead to un- 
physical overscreening and spectral broadening. When 
the QP levels are obtained from the QP equation, the 
Go-dependence is less pronounced because Sgw[Go] is 
evaluated at e° which is consistent with Go- 

The self-consistent GW spectrum is independent of the 
choice of Go, but the number of iterations required to 
reach self-consistency is less when based on Ghf- 



C. Basis set convergence 

In Figs. H] and [5] we show the energy of the three high- 
est occupied molecular orbitals of H2O and CO obtained 
from selfconsistent GW using various sizes of the aug- 
mented Wannier basis. Clearly, the polarization func- 
tions have relatively little influence on the QP energies 
while the first set of additional zeta functions reduce the 
QP energies by up to 0.5 eV. The differences between 
DZP and TZDP are less than 0.15 eV for all the levels 
which justifies the use of DZP basis. 

We have also compared the eigenvalues obtained from 
selfconsistent HF calculations using the DZP augmented 
Wannier basis to accurate HF calculations performed 
with the real-space code GPAW— . Here we obtain a 
MAE of 0.09 eV for the energy of the HOMO level of the 
33 molecules. 



IV. CONCLUSIONS 

As the range of systems to which the GW method is 
being applied continues to expand it becomes important 
to establish its performance for other systems than the 
solids. In this work we have discussed benchmark GW 
calculations for molecular systems. 

The GW calculations were performed using a novel 
scheme based on the PAW method and a basis set consist- 
ing of Wannier functions augmented by numerical atomic 
orbitals. We found that a basis corresponding to double- 
zeta with polarization functions was sufficient to obtain 
GW energies converged to within 0.1 eV (compared to 
triple-zeta with double polarization functions). The GW 
self-energy was calculated on the real frequency axis in- 
cluding its full frequency dependence and off-diagonal 
elements. We thereby avoid all of the commonly used 
approximations, such as the the plasmon pole approx- 
imation, the linearized quasiparticle equation and ana- 
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FIG. 4: (Color online) Convergence of the three highest oc- 
cupied levels of H2O obtained from GW calculations with 
different sizes of the augmented Wannier function basis. SZ 
denotes the Wannier function basis, while e.g. DZDP denotes 
the Wannier basis augmented by one extra radial function per 
valence state and two sets of polarization functions. 



Basis set convergence for CO 



> 

CD 



> 
CD 



> 
CD 



-13.4 

-13.6 

-13.8 

-14 
-15.8 

-16 

-16.2 

-16.4 

-19.4 
-19.6 
-19.8 
-20 



HOMO 


1 1 1 

DZP* B TZP 


- A 

SZP 

1 


/ DZDP^^p J 

/ SZDP/ - 
1 1 1 



HOMO-1 


1 1 1 


T A 




1 


1 1 1 



HOMO-2 


1 1 1 


~ A 

1 


1 1 1 



10 20 30 40 
# Basis functions 



50 



FIG. 5: (Color online) Same as Fig. [4] but for CO. 



lytical continuations from imaginary to real frequencies, 
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and thus obtain a direct and unbiased assessment of the 
GW approximation itself. We found that the inclusion 
of valence-core exchange interactions, as facilitated by 
the PAW method, is important and affect the HF/GW 
HOMO levels by -1.2 eV on average. 

The position of the HOMO for a series of 33 molecules 
was calculated using fully selfconsistent GW, single-shot 
G W , Hartree-Fock, DFT-PBEO, and DFT-PBE. Both 
PBE and PBEO eigenvalues grossly overestimate the 
HOMO energy with a mean absolute error (MAE) with 
respect to the experimental ionization potentials (IP) of 
4.4 and 2.5 eV, respectively. Hartree-Fock underesti- 
mates the HOMO energy but improves the agreement 
with experiments yielding a MAE of 0.8 eV. GW and 
GoWo overcorrects the Hartree-Fock levels slightly lead- 
ing to a small overestimation of the HOMO energy with 
a MAE relative to experiments of 0.4-0.5 eV. This shows 
that although screening is a weak effect in molecular sys- 
tems its inclusion at the GW level improves the electron 
removal energies by 30-50% relative to the unscreened 
Hartree-Fock. The best IPs were obtained from one-shot 
Go Wo calculations starting from the HF Green's func- 
tion where the overscreening is least severe. Very similar 
conclusions were reached by comparing GW, Go Wo and 
HF to exact diagonalization for conjugated molecules de- 
scribed by the semi-empirical PPP model^ 
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Appendix A: The GW self-energy 

Let U denote the rotation matrix that diagonalizes the 
pair orbital overlap Sj^jy = (nij\riki), i.e. U'SU = al. 
The columns of U are truncated to those which have 
corresponding eigenvalues a q < 10 _5 a G 7 3 . We then only 
calculate the reduced number of Coulomb elements 

v w' = ( n <i\y^nM> ( A1 ) 

where n q (r) are the optimized pair orbitals 

n q {r) = Y^ n iA r ) u v,ql V°qi ( A2 ) 

ij 

which are mutually orthonormal, i.e. (n q \n q >) = 8 qq i. 

Determining the GW self-energy proceeds by calculat- 
ing first the full polarization matrix in the time domain 

P£ kl (t) = 2iG<(t)G>*(t), (A3) 

^hi® = P jtlM (A4) 



The factor 2 appears for spin-paired systems from sum- 
ming over spin indices. This is then downfolded to the 
reduced representation 

p i = E v^ u k, p iki u kwv^- (as) 

ij,kl 

The screened interaction can be determined from the 
lesser and greater polarization matrices, and the static 
interaction V qq > , via the relations 

P r (t) = 9(t)(P>(t)-P<(t)), (A6) 

W r {e) = [I- VP r (e)]- 1 V, (A7) 

W>{e) = W r (e)P>(e)W r \e), (A8) 

W<(e) = W>(e) - W r (e) + W r \e), (A9) 

where all quantities are matrices in the optimized pair 
orbital basis and matrix multiplication is implied. We 
obtain the screened interaction in the original orbital ba- 
sis from 

W Sm (*) * E U ii,*V^wf g , , (A10) 

qq> 

which is an approximation due to the truncation of the 
columns of U. Finally the GW self-energy can be deter- 
mined by 

siw,«(*)=iE G ^(*) w S,ii( < ) ( An ) 

kl 

Sgw(') = 0(t) (£> w (i) - £< w (t)) + (A12) 
The exchange and Hartree potentials are given by 

s,,ii = *E^ G w(* = ) ( A13 ) 

kl 

E Hl « = -2i J2 v ^Gti (t = 0) (A14) 

kl 

The Green functions are given by 

G r (e) =[(e + i V )S - H KS + v xc - A m - ^(e)]- 1 

(A15) 

G<{e) = - /fd(£ - M)[G r (e) - G r (e^} (A16) 
G>(e) =(1 - / FD ( £ - n))[G r (e) - G r (e)^} (A17) 

where /fd (s — A 4 ) is the Fermi-Dirac function and Ai>h = 
Eh[G] — Sh[Gdft] is the difference between the GW 
Hartree potential and the DFT Hartree potential. For 
self-consistent calculations, equation (|A3I) - (IA17|) are it- 
erated untill convergence in G. 
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